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Abstract : 

A very simple piijcedure has been developed to fit the first three moments 
of an actual phase function with a three parameter analytic phase function. The 
exact Legendre Polynomial decomposition of this function is known which makes 
it quite suitable for multiple scattering calculations. The use of this function 
can be expected to yield excellent flux values at all depths within a medium. 
Since it is capable of reproducing the glory it can be used in synthetic spectra 
computations from planetary atmospheres. Accurate asymptotic radiance values 
can also be achieved as long as the single scattering albedo tOoSO-9* 


! 
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I, Introduction 

One of the most widely used analytic phase functions In radiative transfer 
theory is that of Henyey and Greenstein. O) Their function is specified in 

j| 

I 

terms of a single parameter, called the asymmetry parameter, usually denoted by g, 
and has a simple Legendre polynomial decomposition. This function is given by 


♦ (li* s) (l-g^)/[An (l-2gp+g^)3/^] 
where p is the cosine of the scattering angle and 

g = <P> =, J \ pc>(p, g) dpd^) 

(1) 


(2) 

li 

^ ai| 

This function has the following normalization 


j 

9.rr 1 

\ \ g) dpd(^ = 1 . 

i 

Now using the generating function for the Legendre polynomials, 

(3) 

namely 

i 

t 

(l-2gp + = “2 

n=0 

and differentiating with respect to p, it is easy to show that 

go 

(4) 

1 

1 

i 

g) = 1 ’2g" (2n+l) Pn (p). 
Att tf=0 

Defining 

^ p 1 

(5) 


<Pjj(p)> = j 3 Pn(p)«(p, g) dpd<i> 
0 

and using the fact chat 

r‘ 

(6) 


]^Pn(p)Pn,(M) dp =2^ 

(7) 


we get 



<Pn(M)> = 

(8) 



One of the majcr disadvantages of this single term phase function is 
that it can't produce the familiar glory which is almost always prevelant 
in realistic phase functions computed from Hie Theory. This is of major 
importance when one is computing synthetic spectra in planetary atmospheres. 
This fact was pointed out by Whitehill and Hansen^^^ to explain the so called 
inverse phase effect on Venus. To produce the glory phenomena requires the 
use of a superposition of two Henyey-Greentstein phase functions. This was 

Co) 

done by Irvine who used 
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■Kwi Si> S2» a) “ a 'Km. gj^) + (l-a)*(p, ^ 2 '^ 

However, Irvine only fit the flifjt moment of the phase function. It is the 

purpose of this paper to show how the parameters g 2 » and a can be computed 

very simply to fit at least three moments of the phase function thus allowing 

more accurate flux and internal radiance calculations. Hereafter the function 

defined by equation (1) will be referred to as OTHG (one term Henyey-Greenstein) and 

that defined by (9) as TTHG (two term Henyey-Greenstein). 

2. Theory 

2.1 Moment Fitting 

Let| us assume that we are given the first three components of the pro- 
jection of a realistic phase function on a basis of Legendre polynomials and 
label them g, h, and t respectively. This is equivalent to knowing the first 
three moments. From equations (5) and (9) it follows that 


g^, g2» «) = 1^21 + (l-a)g2'^)(2n+l) P„(p) (iq) 


and 

n=o 



<Pj^(p)> = + (l-a)g2^ 

<ii) 

We can now set up 

the following three equations 



agj^ + (l-a)g 2 = 8 

(12) 


ag^^ + (l-a)g2^ “ ^ 

(13) 


agj^^ + (l-o)g 2 ^ = t 

(14) 

These equations can readily be solved yielding 



82 = t-hg-[(hg-t)^-4(h-g2)(tg-h2)]l/2/2(h-g2) 

(15) 



(16) 


a = ( 8 - 82 ) 

(17) 


The minus sign was chosen so that g 2 could become negative. If one 
chooses the plus sign for the square root in equation (15) then one obtains 
It should be noted here that there is nothing In these equations to prevent 
|g 2 l>l. Although this case is the exception rather than the rule we will show 
how to deal with it later. 
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2.2 Internal Radiance from Eigenvalue Expansion 
(A) 

Van de Hulst * derived an expression for the asymptotic internal radiance 
in a homogeneous medium in terms of a perturbation expansion involving an eigen- 
value along with g and h. We will extend the perturbation expansion to also 
include t and thus show that a better fit can be obtained. We will adopt the 
notation of van de Hulst to avoid any confusion to the reader. The basic 
assumption is that the intensity deep inside a medium and far away from a lower 
boundary can be expressed as 

I (Tp) a (P), - 1 < p 1 1 <18) 

where t is the optical depth measured from the top of the medium and k is the 
smallest eigenvalue of the equation 


where 


(1-kp) P(p) =1^ J h(p, p') P(p'') dp" 

I 

J-3LTI 

h(p, p”") = J<!'[pp^ + (l-p2)^ (X-p"2)% cos(i}>-<|i")l d(|) 


(19) 


( 20 ) 


and Wq is the single scattering albedo. 

(5) 

This assumption hap been checked numerically by Plass et al. and was 
shown to be quite valid for oIq 1. The exponential behavior has also been 
demonstrated exactly by a solution to the one dimensional problem (Kattawar 
and Plass ^^^) . 

Van de Hulst obtains the following series expansions for Wq and P(p); 


10 = 1 - 
"o 


3(l-g)^ 45(l-g)^(i-h)^ + ; 


( 21 ) 


P(„) . 1 + + Oft') 


( 22 ) 
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Since we have demons crated a method of fitting a phase function up to the 
term for we must extend the perturbation series for P(y) up to this 

term. This can be done in a rather straightforward manner and we obtain 


P(y) - 1 + + 5(1-8) 

+ OOc**) (23) 

The coefficient of in will involve <P^(y)> and we will not include it. 

It is thus clear that by using the three parameter pha^iS function that more 

i 

accurate’ internal radiances will be obtained since the OTHG only fits g. It 

I 

should be noted, however, that all higher values of <Pj^(p)> are given by g 
for the OTHG. In order to carry out the expansion given by equation (23) we 
must be able to determine the eigenvalue k. 

2.3 Eigenvalue Computation 

The eigenvalues can be computed from equation (19) . The proof that the 

(7) 

eigenvalues are real was given by Kuscer and Vidav . The integral was dis- 
cretized by using a twenty eight point Lobatto quadrature. Now equation (19) 
can easily be manipulated so that one can solve for k^ rather than k since it is 
well known that the eigenvalues occu» in pairs (±k) . The program was tested by 
using the isotropic phase function (g-o) and computing the eigenvalues for 

( 8 ) 

several values of The results were compared to those given by van de Hulst 

for this case and agreement was achieved to every figure given. 

3. Numerical Results 

The first thing to be tested was how well the eigenvalues computed by the 
OTHG and the TTHG compared to the eigenvalues of an actual phase function. The 

. (9) 

phase function was computed using the Haza L distribution of Deirmendjian 
which is a modified gamma distribution of the form 


n(r) 


2 -15. 11 86 rH 

are 


(24) 


The refractive index used was N = 1.55-O.Oi typical of atmospheric aerosols for 
a wavelength of 0.55 pm. The resulting phase function la shown in Figure 1 
along with the OTHG and TTHG which were computed from the appropriate values of 
<P^>, The values of (n ® li 2, 3) were computed in two ways. First a 

numerical integration was performed over the phase function computed at very 
small angular Intervals. Secondly these values can also be obtained directly 


from the A coefficients (See Kattawar et al. 

(n+1) (n+2) 


( 10 ) 


) by 


<P > 
n 


4tt I (n+1) (n+2) .1 , 6n^ + 6n-4 1 

2n+l (2n+3) (2n+5) ^ n+2 (2n-l) (2n+3) 


n 


n(n-l) 


+ 2(n-M) , 3 . 

(2n-3)(2n-l)'^ n-2 (2n+3) ' n+1 (n-1) 


2n 3 T 
-1) ^ n-lj 


(25) 


Excellent agreement was obtained between the tx^o methods. It should be noted 
from Figure 1 that the TTHG not only provides a better approximation to the 
glory (0“18O°) but also a better approximation to the diffraction peak (0=0°), 

Table 1 gives a comparison between the eigenvalues computed from a realistic 
phase function and those computed using a TTHG and a OTHG approximation to it. 

It should be emphasized, that only the smallest eigenvalue was chosen since this 
is the one that dominates the solution for large t. It should be noted that the 
TTHG provides excellent eigenvalue agreement for all values of o)q used. 

As a further test we wanted to consider the efficacy of the method in de- 
termining internal radiance in optically deep homogeneous media. Our first test 
employed the well known Rayleigh phase function. Since g = 0 for this function, 
the OTHG reduces to isotropic scattering; however, h = 0,1 and the TTHG produces 
the curve shown in Figure 2, This is a much more realistic description of Rayleigh 
scattering than the OTHG. The asymptotic internal radiance for Rayleigh scat- 
tering has already been calculated by Plass, et al.^^^ The eigenvalue for the 
case where Wq = 0.9 was computed both from equation (19) and by fitting the actual 
exponential fall off. Both calculations gave the same value for the eigenvalue 
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namely k = 0.5232. The values of k obtained by using the OTHG and the TTHG 
were 0.525A and C.5232, respectively. A comparison of the upward and downward 
internal radiance is shown in Figure 3. The radiance values presented for the 
OTHG and the TTHG were computed from equation (23) using the appropriate eigen- 
values. As can be seen from the figure, the OTHG provides a better fit for oiy<0.2 
for the upward radiance but the TTHG is better for 0.2<u£l.0 For the downward 
radiOTce the TTHG is better everywhere. It should also be noted that for values of 
0 ) >0.9 one can obtain quite excellent asympotic internal radiances without doing very 

° I 

large scale computer calculations by simply using equation (23) with appropriate values 
of k, g, h, and t. 

Ve next computed a phase function using the size distribution in equation (24) and a 

wavelength of 0,55 pm for a refractive index of N^l.SS-O.OOli, This curve along with 

the OTHG and TTHG are shown in Figure 4. The asymptotic normalized radiances -are shown 

in Figure 5. The difference between the actual asymptotic radiance and that computed 

with the TTHG using the parameters k, g, h, t, in equation (23) is almost imperceptible 

on the scale of the figure. The OTHG internal radiances computed with the parameters k, 

2 3 

Si g f 8 is not quite as good, A more dramatic comparison can be seen in the following 
exAmplev The size distribution was held fixed and the refractive index was changed to 
N“l. 55-0. 01271. This gave the following set of curves shown in Figure 6. In Figure 7 we 
show^a comparison of the asymptotic radiances computed using the matrix operator theory 
with the actual phase function compared to that for the TTHG and OTHG computed 
using equation (23). For the TTHG the parameters k, g, h, and t were used which 
agree with the actual phase function and for the OTHG the corresponding parameters 
were k, g, g^ , and g^. It is easily seen that the TTHG provides a better approxi- 
mation to the upward radiance. There is always more disagreement in the upward 
radiance as compared to the downward radiance, since the upward is derived from 
the downward. Although the curves shown for the TTHG and the OTHG were computed up 
to order k * the expansion can be extended to higher orders until convergence is 
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reached. This would give the true asymptotic limit since all values of <P^> 
are known for these t\Jo functions. 

To test the TTHG on a phase function indicative of terrestrial clouds, we 
xised the size distribution 


. , 6 -1.5r 

n(r) are 


(26) 


with N » 1.33-O.Oi for X = 0.77 pm. The following values were obtained; 
g - 0.8479, h “ 0.7763, and t = 0.6523. Fitting these three values the TTHG 
gives g. = 0.8792, g„- 0.9835, and a » 0.9832. These results are presented in 

Figure 8 along with the OTHG. It is easily seen that the TTHG greatly over- 

estimates the glory peak(e = 180®) even though this function has the same values 
of g, h, and t as the actual phase function. If one, however, wants to tolerate 
a small error in t and fit the glory peak this can be done in the following way. 

We will still require that the TTHG fit both g and h, however, we will relax the 

requirement on fitting t in favor of fitting a particular value for the glory. 
Therefore, 


Eliminating g^ we get 


+ (l’-a)S2 " ® 

(27) 

ag^ + (l-a)g2 = h 

(28) 

^ G - [o(H-G^)/(l-a)]^ 

(29) 

= [G-(l-a)g2]/a 

(30) 


Therefore by adjusting the ratio a/(l~a) we can reproduce the glory peak at the 
expense of introducing an error in t. This is shown in Figure 8, The exact 
value of t is 0.6523 and for the TTHG which fits the glory peak it is 0.6659 
which is a 2 percent error. This stculd be contrasted to the OTHG wViich has a 7 
percent error in h and a 6^2 percent error in t. The major difference in the two 
TTHG phase functions only occurs for scattering angles greater than 110®. This 
demonstrates that the third moment of a phase function is quite sensitive to 
shape and one must be very careful in using phase functions which appear to 
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"look” similar to real phase functions. The technique just described is what 
one would use when a fit to g, h, and t produces a value of Of'® would 

then have to sacrifice a fit to t in favor of a fit to the glory. 

In order to study the effects of the approximate phase functions on the 

reflected intensity, we chose the following model; 

, X 6 -12r (31) 

n(r) are ' 

with a refractive index of N = 2,25-0,0061 for a wavelength of 0,565 vin. This 
phase function along with the OTHG and TTHG are shown in Figure 9. The single 
scattering albedo » 0,8943, In order that there be no scaling in optical 
depth we considered a homogeneous semi-infinite atmosphere. The matrix operator 
^program was used to calculate the reflected radiances for the three phase functions. 
In Figure 10 we sec. a comparison of the reflected radiance for = 1,0 as a 

function of the cosine of the zenith angle. What is particularly striking is the 
degree to which single scattering features are carried over into the multiple 
scattering results. For example from Figure 9 we see chat the TTHG crosses the 
actual phase function curve for angles > 90® at two points corresponding to p = 

0,24 and 0,91. These correspond very closely to the crossover points from Figure 
10 for Che multiple scattering results. The values in parentheses are the irradi- 
ance values and it can be seen that the TTHG gives only a 1 percent error in the 
irradiance while the OTHG gives a 7 percent error. Although Che TTHG produces 
good values for the irradiance, the reflected radiance shows large discrepancies. 
Therefore if one needs precise values of the radiance then use of the analytic 
phase function can give errors which are quite large. Figure 11 is the same case 
as Figure 10 except = 0,1182. The radiance in the principal plane is shown 
(({i = 0° and 180®). Again large discrepancies exist bem^een the actual phase 
function results and the two analytic phase functions, however the flux value for 
the TTHG is in error by only 1 percent while the OTHG is in error by 5.5 percent. 
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4. Conclusions 

The use o£ this function can be expected to yield quite good asymptotic 
radiances in a medium whose single scattering albedo Is > 0.9. We have shown 
that this function can yield quite accurate eigenvalues which govern the rate of 
decrease of the radiance and irradiance at asymptotic depths. This also Implies 
its effectiveness in the computation of heating rates. It can also produce the 
familiar glory so often seen in realistic phase functions. This should be quite 
useful in the computations of synthetic spectra from planetary atmospheres where 
a parameterization of the phase function can be Judiciously used. Although the 
values of the reflected radiance con be In error by a large percent, the reflected 
irradiance values can be quite accurate. Therefore this function can be effect- 
ively used whenever radiance values averaged over certain intervals. This 
Is precisely the case when we are computing synthetic spectra from various 
portions of a planetary atmosphere where averages are performed over portions 
of a planet covered by the spectrograph slit. 
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Figure Captions 


Fig. 1 


Fig. 2 


Fig. 3 


Fig. /, 
Fig. 5 

Fig. 6 
Fig. 7 
Fig. 8 

Fig. 9 

Fig. 10 


Fig. 11 


Phase function as a function of scattering angle for a 
size distribution n(r) a r e" * with refractive 

Index N - 1.55 •* 0.01 at a wavelength of X"0.55pm. TTHG 
Is the two term Hiinjcy -Greens tein analytic phase function 
obtained by fitting the first three moments of the actual 
phase function and OTHG Is the one term Henyey Greens teln 
phase function obtained by fitting only the first moment. 

Some as Fig. 1, except the actual phase function Is that 
for Rayleigh scattering. 

Upward and downward normalized asymptotic Internal radi- 
ances, The values of the single scattering albedo u^were 
computed from equation (21) using the eigenvalues k com- 
puted from equation (19). These results are for the phase 
function used in Fig, 2. 

Same as Fig, 1 except N “ 1.55 - 0,0011 

Same as Fig. 3 except phase function used was that of 

Fig. 4. 

Sames as Fig. 1 except N = 1,55 - 0.01271 

Same as Fig. 3 except phase function was that of Fig. 6. 

^ ^ 6 -1.5r with N « 1.33 - 0.01 

Same as Fig, 1 except n;r) are 

for X - 0.70pm. 

6 — 

Same as Fig 7 except n(r) a r e with N = 2,25 - 0.006i 
for X ■*= 0.565pm. 

Reflected radiances as a function of the cosine of the zenith 
angle for a homogeneous semi-infinite atmosphere scattering 
according to the phase function for Fig. 9. Cosine of 
incident angle <=■ 1.0. 

Same as Fig. 10 except = 0.1882 and the results are for the 
the principal plane <{» = 0“ and 180®. 
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TABLE I 


Comparison of eigenvalues for a real phase function 
and those computed from the OTHG and the TTJIG. 
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'‘tthg 

% Error 
for TTHG 

^OTHG 

% Error 
for OTHG 

0.9999 

0.01015 

0,01015 

0 

0,01015 

0 

0.999 

0.03209 

0.03209 

0 

0.03209 

0 

0.99 

0.1016 

0.1016 

0 

0.1017 

0.1 

0.95 

0.2284 

0,2284 

0 

0.2296 

0.23 

0.90 

0.3251 

0.3251 

0 

0.3281 

0.92 

0.85 

0.4006 

1 0.4006 

0 

0.4054 

1.20 

0.80 

0.4652 

0.4652 

0 

0.4718 

1.42 

0.75 

0.5227 

0.5227 

0 

0.5310 

1.59 

0.65 

0,6237 

0.6238 

0.02 

0.6347 

1.76 

0.55 

0.7115 

0.7116 

0.01 

0.7241 

1 

1.77 

0.45 

0.7890 

0.7892 

0.02 

0.8022 ' 

1.67 

0.35 

0.857S' 

0.8578 

0.03 

0.8701 

1.47 

0.25 

0.9166 

0.9173 

0.08 

0.9275 

1.19 

0.15 

0.9648 

0.9657 

0,09 

0.9724 

0.79 

0,05 

0.9958 

0.9965 

0.07 

0.9978 

0.2 
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Figure 10 











